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1 . 0   Introduction 

The  problem  of  interpolating  or  approximating  data  from 
scattered  measurements  arises  in  many  areas  of  science  and 
engineering.  The  importance  of  the  problem  has  resulted  in  a 
large  number  of  methods  for  solution  of  the  problem,  as  has  been 
noted  in  surveys  by  Schumaker  [55],  Barnhill  [4],  and  more  re- 
cently by  Franke  [19].  The  problem  can  be  described  very  easily. 
Given  data  ( x . , y . ,  f .)  ,  i=l,...,N,  construct  a  (smooth)  function, 
F(x,y)  such  that  F(x.,y.)  -  f . ,  i=l,...,N.  No  assumption  is  made 
regarding  the  spacing  of  the  independent  variable  data.  More 
generally,  especially  when  the  data  are  subject  to  errors,  one 
may  wish  to  relax  the  interpolation  condition  and  approximate  the 
surface.  The  problem  has  an  obvious  generalization  to  more 
independent  variables.  The  existence  of  many  methods  for  such  a 
surface  is  due  to  the  many  sources  of  data  and  the  great  number 
of  possible  dispositions  of  data  points.  Some  advice  on  how  to 
proceed  for  given  data  is  contained  in  Sabin  [50].  It  is  clear 
that  no  one  method  is  satisfactory  for  all  cases. 

In  the  last  few  years  a  number  of  advances  have  been  made. 
Because  of  space  and  time  limitations,  this  paper  will  be  con- 
fined to  a  few  of  those  which  I  feel  to  be  particularly  important 
or  interesting.  The  first  area  concerns  the  mathematical  under- 
pinnings of  the  multiquadric  method  of  Hardy  [22] .  This  method 
has  previously  been  noted  to  work  well  in  a  variety  of  cases  (see 
Franke  [19]  and  Kansa  [29],  for  example),  but  until  recently 
little  was  known  in  terms  of  a  mathematical  theory  for  the 
method.    These  developments  will  be  reviewed  in  section  2.    The 


construction  of  surfaces  with  tension  parameters  or  satisfying 
constraints  seems  to  be  desirable  based  on  the  importance  of  such 
ideas  in  the  univariate  case,  and  recent  results  are  discussed  in 
section  3.  Finally,  since  the  amount  of  data  is  sometimes  far  in 
excess  of  what  is  required  to  define  a  sufficiently  accurate 
surface,  the  problem  of  surface  approximation  and  selection  of 
subsets  on  which  to  base  interpolating  (or  approximating)  sur- 
faces has  become  important,  and  this  topic  will  be  treated  in 
section  4 . 

There  have  been  other  developments  which  are  important,  but 
will  not  be  treated  here.  These  include  convergence  of  Shepard's 
method  (see  Farwig  [17]),  interpolation  on  the  sphere  (see  Lawson 
[32],  Wahba  [63],  Renka  [49],  Ramaraj  [48],  and  Nielson  and 
Ramaraj  [46] ) ,  and  some  new  implementations  of  "finite  element" 
methods  and  related  things  such  as  derivative  estimation  and  new, 
higher  order  elements  (see  Alfeld  [2-3],  LeMehaute  [33],  [34], 
Sablonniere  [51]).  Certain  methods  based  on  statistical  ideas 
(especially  Kriging  in  geology  and  Optimum  Interpolation  in 
meteorology)  continue  to  be  the  focus  of  much  effort  in  those 
particular  disciplines.  Relevant  recent  references  include 
Journel  [28]  and  Thiebaux  [59]. 

2.0   Multiquadrics  and  Related  Methods 

The  multiquadric  (MQ)  method  was  proposed  by  Hardy  [22] ,  and 
he  has  investigated  the  fitting  properties  of  the  method  when 
applied  to  data  from  various  sources  in  a  series  of  papers  that 
have  appeared  since  that  time  [23-25].  The  scheme  is  quite 
simple   to   implement   and  is  reasonably  efficient   in   terms   of 


computer  resources  provided  the  number  of  points  is  not  large.   A 

2     2  1/2 
basis   function  (a  quadric)  B.  =  (d.  +  r  )  '       is  associated   with 

V  J 

J.  1» 

the  j    data  point.   Here  d.  is  the  distance  from  the  point  (x,y) 

to  the  j    data  point  (x.,y.)  and  r  is  a  parameter  in  the  method. 

Note   that  each  of  the  basis  functions  is  a  radial  function   with 

respect   to  the  data  point,   and  that  they  are  all  translates   of 

each  other,  which  overcomes  the  usual  problem  of  running  afoul  of 

the   Haar   theorem   regarding   interpolation   in   more   than   one 

variable.   Now,  a  linear  combination  of  the  functions 

N 
F(x,y)  =  y\..  B.(x,y) 

3  =  1 
is   required   to  interpolate  the  given   data.    This   yields   the 

system  of  equations, 
N 


E 


3j(xi,yi)aJ  =  f.   ,  i=l,2....,H 


The  existence  of  the  interpolant  is  dependent  on  the  nonsingu- 
larity  of  the  matrix  {B  .(x  .  ,  y  .  )  }  .   Other  authors  have  also  inves- 

J  -*-  -L 

tigated  the  scheme,  generally  from  an  empirical  point  of  view 
[19],  [27],  [57].  These  authors  have  found  that  the  method  is 
quite  adept  at  yielding  good  approximations,  in  many  situations 
better  than  any  other  method.  In  addition  to  this  method,  the 
use  of  the  reciprocal  of  the  above  basis  function  leads  to  the 
"reciprocal"  MQ  method. 

The  MQ  method  is  one  of  a  class  that  I  have  previously 
called  "global  basis  function  methods",  and  which  others  have 
called  "kernel"  methods.  In  general  the  approximation  takes  the 
form 


N  M 

F(x,y)  =  J\.   B.(x,y)  +  £V  p.(x,y)  , 
j=l  j  =  l 


where  {p.}  is  a  set  of  monomials  of  degree  <m.   The  equations  for 

*J 

the  coefficients  in  such  methods  may  be  written  in  the  form 
N  M 

I]Bj(Xi'7i)aJ  +  SPi(Xi'7i)bi  =  fi'  i=1'  ■■■'*   • 
j  =  l  .3  =  1 

N 

SPi(Xj'yj)aj    =  °    '  t=1'  •••   'M   • 
j  =  l 

The  first  of  the  equations  require  interpolation  to  the  data  by  a 
linear  combination  of  the  basis  functions  B.(x.y)  plus  the  M 
polynomial  terms,  P«(x,y),  while  the  last  set  requires  the  coef- 
ficients  to  satisfy  a  certain  constraint,  which  as  we  shall  see, 
may  related  to  the  conditional  positive  definiteness  of 
(B.(x.,y.)},  and  serves  to  guarantee  exactness  for  the  set  of 
polynomials  {p-(x,y)}.   In  matrix  form,  we  may  write 

\E'  0/\b/    W  . 


2.1   The  Multiquadric  Method:   Theory 

The  intriguing  aspect  of  the  method  is  that  until  recently, 
very  little  was  known  in  terms  of  a  mathematical  basis  for  the 
efficacy  of  the  method,  even  whether  or  not  the  coefficient 
matrix  was  possibly  singular.  As  recently  as  1983  at  The  Inter- 
national Symposium  on  Surface  Approximation,  in  Gargnano ,  Italy, 
I  proposed  as  a  conjecture  the  inequality 

(-l)N"1det{Bj(xi>y.)}  >  0  . 
As  it  turned  out,  Charles  Micchelli  promptly  heard  of  the  conjec- 
ture,  and   subsequently  proved  it,   and  along  the  way,   theorems 


which  answered  some  other  questions  as  well. 

To  discuss  the  results  of  Micchelli  [41]  it  is  necessary  to 
define  some  terms  and  give  some  background  information.  Because 
of  the  generalization  to  s-dimensional  space,  for  this  discus- 
sion,  points  in  R  will  be  denoted  as  (possibly  superscripted) 
boldface  vectors  rather  than  subscripted  coordinates. 

Definition:   A  continuous  function  F(t),  defined  on  [0,oc)  is  said 

to  be  conditionally  (strictly)  positive  definite  of  order  k  on  R 

if  for  any  distinct  points  x  ,  x  ,  .  .  .  ,  x   in  R  ,  and  scalars  c.  , 

c0,  . . . ,  c   such  that 
Z        n 

n 


^Tc.ptx1)  =  0 


i=l 

for   all  polynomials  p  over  R   of  degree  <k,   the  quadratic   form 


n    n 


ZEcicJF(,ixi  - 


xj||2) 


i=l  j=l 

is  (positive)  nonnegative. 

Let  the  class  of  functions  which  are  conditionally  positive 
definite  of  order  k  over  R  be  denoted  by  P^(R  ) .  and  the  class 
of  functions  which  are  conditionally  positive  definite  of  order  k 
over  Rs ,  for  all  s,  by  P,  .  Further,  recall  that  a  function  F  is 
completely  monotonic  on  (0,°o)  if  it  is  in  C  (0,°o)  and 
(-l)mF(m) (t)  >  0  for  all  t>0  and  m=0,l,2,...  . 

Theorem  1:   F  is  in  Pi  whenever  F  is  continuous  on  [0,oo)  and 

k  ( k) 
(-1)  Fv   (t)  is  completely  monotonic  on  (0,oo). 

This  theorem  is  due  to  Schoenberg  [54]  for  order  k=0 .    Fur- 
ther, if  F  is  not  a  polynomial  of  degree  <k,  and  if  the  points  in 


"the   definition   of  completely  monotonic  are  distinct,   "then   the 

2    1  /2 

quadratic   form  is  positive.    Take  F(t)  =  ( r  +t )    ,   and    note 

that   for  m>l,   F(m)(t)=C  ( -1 )m~1 ( r2+t ) ( 1_2m)/2    where  C    is   a 

m  m 

positive   constant.    This  special  case  leads  to  the   coefficient 

matrix   for   the   MQ  method,   which  is  seen  to   be   conditionally 

positive  definite  of  order  k=l. 

When  the  constant  r  is  taken  "to  be  zero,  each  basis  function 

is  the  upper  half  of  a  cone,   hence  the  interpolant  is  not  smooth 

at   the  data  points.    This  particular  interpolant  is  almost   the 

"multiconic"  method  of  Duchon  [13].    The  difference  is  that   the 

multiconic   approximation   is   consistent   with   the   concept   of 

conditional   positive   definiteness  of  order   one.    The   overall 

approximation  takes  the  form  of  a  linear  combination  of  the  basis 

functions   plus  a  constant,   with  the  additional  constraint   that 

the  coefficients  satisfy, 
n 

Eao  =  ° 


This  constraint  is  easily  seen  to  guarantee  precision  for  con- 
stant functions.  This  special  case,  along  with  Theorem  1,  is 
convincing  evidence  that  a  constant  and  the  corresponding  con- 
straint should  be  included  in  the  MQ  approximation.  In  a  practi- 
cal vein,  however,  my  own  limited  tests  have  indicated  that  the 
accuracy  of  the  method  is  not  always  helped,  and  may  be  hindered, 
by  including  the  constant. 

The  "reciprocal"  MQ  method  uses  2F'(t)  for  the  basis 
function,  and  thus  the  theorem  shows  that  the  coefficient  matrix 
for  that  scheme  is  positive  definite. 
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Theorem  2:  Let  l=[s/2]-k+2  be  a  positive  integer.  Then  for  any 
function  defined  on  (0,°o  )  such  that  ( -1 )  % ^  J ' ( t )  is  nonnega- 
tive,  nonincreasing,  and  convex  for  j=0,  1,  . . . ,  1-2  on  (0,qo)  (if 
1=1,  we  require  only  that  it  be  nonnegative  and  nonincreasing), 
it  follows  that  F(\/t)  is  in  Pk(Rs). 

Theorem   3:    Suppose  F'  is  completely  monotonia  but  not  constant 

on  ( 0 ,  oo  )  ,   F  is  continuous  on  [0,^)  and  positive  on  (0,°°).   Then 
for  any  distinct  vectors  x,x,  ...,x   in  R   (s  arbitrary) 
(-l)n_1det{F( ! |xi-xj | |2)}  >  0  . 

This  last  theorem  proves  the  conjecture  about  the  coeffi- 
cient matrix  in  the  MQ  method.  However,  the  theorem  also  yields 
similar  results  for  interpolation  by  other  sets  of  radial  basis 
functions .  Some  of  those  methods  are  known  by  other  results  to 
always  lead  to  nonsingular  systems  of  equations,  such  as  the 
"thin  plate  splines"  of  Duchon  [12-13]  (see  also  Harder  and 
Desmarais  [21]  and  Meinguet  [37-40]),  which  are  known  to  be 
nonsingular  because  of  their  semi-Hilbert  space  setting.  Others, 
such  as  the  basis  function  log( 1+ M x-xJ I  i  )  suggested  by  Dyn  [16] 
are  also  seen  to  be  positive  definite  of  order  k=l. 

Another  result  of  interest  in  specific  applications  has  been 
obtained  by  Hardy  and  Nelson  [26].  This  result  shows  that  the  MQ 
method  has  a  basis  for  approximation  of  geodetic  and  gravita- 
tional anomalies.  The  connection  comes  through  the  representa- 
tion of  the  disturbing  potential  at  point  i  due  to  anomalous 
disturbances  as  a  three  dimensional  integral, 


J\7 


P  dV  , 

V 


where  d.  is  the  distance  from  the  point  i,  p  is  the  Laplacian  of 
the  disturbing  potential,  and  V  is  the  volume  of  interest.  What 
Hardy  and  Nelson  showed  is  that  the  MQ  method  can  be  interpreted 
as  a  quadrature  approximation  to  the  above  integral,  one  that  is 
required  to  yield  the  correct  result  at  certain  (the  data) 
points.  While  this  particular  result  may  say  little  about  the 
scheme  as  a  general  interpolation  scheme,  it  is  interesting  that 
such  an  interpretation  exists  for  one  of  the  early  the  uses  of 
the  method. 

The  role  of  the  parameter  r  in  the  MQ  method  has  long  puz- 
zled investigators.  The  original  interpretation  given  by  Hardy 
was  a  three  dimensional  one:  His  applications  were  actually  in 
3-space,  and  this  value  simply  represented  the  z  coordinate  of 
the  locations  of  the  disturbing  (point)  potential.  Given  that 
the  method  performs  so  well,  it  seems  likely  that  the  MQ  approxi- 
mation can  be  described  as  approximation  in  some  Sobolev-like 
space  similar  to  that  for  the  multiconic  method.  A  recent  devel- 
opment by  Madych  and  Nelson  [35]  gives  this  result.  The  MQ 
method  does  minimize  a  certain  pseudonorm  involving  a  weighted 
integral.  Further,  the  result  applies  to  other  similar  schemes, 
and  thus  may  also  lead  to  ways  of  deriving  other  approximations 
with  desired  properties  through  explicit  minimization  of  weighted 
pseudonorms .  If  practical  for  computational  purposes,  such  re- 
sults could  have  far  reaching  implications  in  applied  scattered 
data  approximation. 

2.2  Multiquadric  Method:   Practice 

In  addition  to  the  developments  on  the  mathematical   aspects 


of  the  MQ  method,  some  progress  has  been  made  in  attempting  to 
solve  such  systems  of  equations  by  iterative  means.  A  series  of 
investigations  by  Dyn  and  coworkers  [15-16]  have  studied  the 
problems  of  fitting  scattered  data  with  linear  combinations  of 
radial  basis  functions.  The  basic  idea  comes  from  the  fact  that 
thin  plate  splines  have  basis  functions  which  are  the  fundamental 
solution  of  the  biharmonic  equation.  The  MQ  basis  functions  are 
solutions  of  Laplace's  equation  (in  three  dimensions).  Thus, 
when  certain  finite  difference  approximations  to  the  iterated 
Laplacian  are  applied  to  the  equations,  the  resulting  equations 
tend  to  have  a  large  diagonal  term,  which  then  yields  a  system  of 
equations  amenable  to  solution  by  iterative  schemes.  This  pro- 
cess may  be  thought  of  as  applying  a  conditioning  operator  to  the 
system  of  equations. 

The  idea  of  the  conditioning  operator  is  to  transform  that 
part  of  the  system  of  equations  involving  A  into  an  equivalent 
one  which  is  better  conditioned,  perhaps  even  diagonally 
dominant.  In  addition  to  transforming  A  suitably,  the  operator 
is  constructed  to  annihilate  E.  This  yields  a  singular  system 
CAa  =  Cf ,  which  has  a  unique  solution,  subject  to  the  constraint 
E'a  =  0. 

Construction  of  the  conditioning  operator  involves  triangu- 
lation  of  the  convex  hull  of  the  data  point  locations.  Finite 
difference  operators  are  then  derived  which  approximate  the  nec- 
essary derivatives  on  the  basis  of  the  function  behavior  at  the 
vertices  of  the  triangles.  Some  care  is  necessary  to  ensure  that 
the  operators  annihilate  E,  which  corresponds  to  annihilating  all 
polynomials  of  degree  <m  on  the  given  set  of  data  points. 


The  results  of  applying  these  ideas  to  only  a  limited  number 
of  irregular,  but  "quasi-regular"  grids  are  reported  in  [16]. 
Nonetheless,  the  results  are  very  encouraging,  with  the  condition 
number  of  the  matrix  being  decreased  by  factors  of  up  to  200  or 
more  for  the  MQ  method  with  up  to  121  data  points.  In  addition 
to  the  MQ  method,  the  ideas  are  applicable  to  approximation  by 
thin  plate  splines  and  other  radial  functions,  and  these  are 
reported  on  as  well. 

One  further  interesting  aspect  of  iterative  schemes  based  on 

this   conditioning  is  that  for  the  MQ  method,   shifted   logarithm 

2      2 
(log(d    +  r  )  )  basis  functions,   and  shifted  thin  plate   spline 

2       2       2     2 
((d    +   r  )log(d   +  r  )  )  basis  functions  certain  nice   spectral 

properties  of  the  matrix  CA  seem  to  occur.  In  particular,  it  was 
noted  that  computationally,  the  "rough"  eigenvectors  correspond 
to  the  smaller  eigenvalues.  Certain  iterative  schemes  will  re- 
move those  components  quickly,  leaving  the  surface  corresponding 
to  the  approximate  solution  (before  iteration  to  convergence)  as 
a  smooth  one.  Thus,  terminating  the  iterative  scheme  prematurely 
corresponds  to  a  smoothing  scheme .  Numerical  evidence  and  exam- 
ples are  given  by  Dyn,  Levin,  and  Rippa  [16].  See  [14]  in  this 
Proceedings  for  more  recent  results. 

3.0   Surfaces  with  Tension  or  Constraints 

In  the  univariate  case,  development  of  curves  with  tension 
parameters,  constrained  approximation,  and  monotone  and  convex 
approximation  is  at  a  high  level.  In  each  situation,  several 
reasonable  schemes  for  obtaining  such  approximations  are  readily 
available.    In  two  or  more  independent  variables,   for  scattered 
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data,  the  situation  is  not  nearly  so  well  developed.  Indeed,  the 
idea  of  exactly  what  characterizes  monotone  behavior  of  scattered 
data  has  not  yet  been  clearly  given.  Nonetheless,  there  have 
been  some  noteworthy  developments  in  the  general  area  of  surfaces 
with  tension  and  surfaces  satisfying  certain  constraint  relation- 
ships, both  from  theoretical  and  computational  points  of  view. 

3.1   Tension 

A  generalization  of  splines  under  tension  to  the  scattered 
data  case  was  given  by  Nielson  and  Franke  [45].  The  basic  idea 
is  that  of  approximation  of  the  the  surface  by  a  "finite  element" 
method,  such  as  proposed  by  numerous  authors,  among  them  Dooley 
[10],  Akima  [1],  Lawson  [31],  Nielson  [43],  and  LeMehaute  [33]. 
The  first  step  of  this  process  consists  of  trianguiation  of  the 
x-y  data  in  some  reasonable  way  (e.g. ,  using  the  max-min  angle 
criterion).  A  certain  "finite  element"  function  is  assumed  over 
each  triangle.  Ordinarily  one  wants  a  smooth  surface,  so  at 
least  C  functions  are  usually  used.  Depending  on  the  element 
chosen,  certain  derivatives  must  be  estimated  from  the  data.  It 
is  this  stage  of  the  process  which  is  crucial  to  the  effective- 
ness of  the  method  (see  Nielson  and  Franke  [44]). 

The   incorporation   of   tension  into   the   approximation   is 

achieved   in   the  following  manner.    Let  the  vertices  (x-y   data 

pairs)  of  the  trianguiation  be  denoted  by  the  set  {V. },   an  index 

set  of  triangle  edges  joining  V.  and  V.  by  N  ,   and  the  edges   of 

i       J      e 

the   trianguiation  by  the  set  {e . .  :ij£N  }.    The  ordering  of   the 
edges  relative  to  the  ordering  of  the  data  points  is  unimportant, 
but   each   edge  must  appear  only  once  in  the  set.    Let  E  be   the 
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union  of  the  set,  of  edges.    Define  the  set,  of  functions  C(E)   as 
those   which  are  C   over  the  convex  hull  of  the  {V.},   restricted 
to  E.    Let  the  vector  <a.   .>  be  given,   with  nonnegative   compo- 
nents ,   each   corresponding  to  a  tension  parameter  for   the   ij 
edge.   Now  define  the  pseudonorm 


8-(»)  =  Tl       /"r(f-f)2  ♦  a?.(fE)2]  de.  . 

ij*N   -'e-. 

e   ij 


There  is  a  unique  minimizer  of  S^(F)  over  all  curve  networks  in 
C(E)  which  interpolate  the  data.  That  minimizer  is  the  function 
which  is  a  piecewise  Hermite  exponential  (a  linear  combination  of 
1,  u,  axp(Q.  .u),  exp(-a. .u) ,  where  u  represents  Euclidean  dist- 
ance  along  the  edge,  that  takes  on  prescribed  value  and  deriva- 
tive conditions  at  the  endpoints ) ,  and  satisfying  a  certain 
system  of  (sparse)  linear  equations  for  the  partial  derivatives 
at  the  vertices.  Asymptotic  properties  are  as  anticipated.  For 
simplicity,  take  all  tension  parameters  to  have  the  same  value, 
and  then  as  tension  goes  to  zero,  the  curve  network  approaches 
the  curve  network  minimizing  the  corresponding  functional,  as 
developed  previously  by  Nielson  [43].  As  the  tension  becomes 
large,  the  curve  network  approaches  the  piecewise  linear  network 
over  the  edges . 

After  obtaining  the  curve  network  over  the  edges  of  the 
triangulation ,  the  surface  is  then  completed  using  a  C  blending 
method  on  the  individual  triangles.  It  is  desirable  to  propogate 
the  effects  of  tension  into  the  interior  of  the  triangles,  so  the 
previously  known  cubic  blending  techniques  are  inadequate.  The 
method   used   was   an  extension  of  the   side-vertex   scheme   (see 
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Nielson  [42])  where  the  radial  projectors  were  taken  to  be  expo- 
nential Hermite  functions.  The  net  effect  of  the  tension  becom- 
ing large  is  that  the  surface  tends  toward  the  the  surface  which 
is  piecewise  linear  over  the  triangulation. 

The  limiting  behavior  of  the  above  construction  is  not  the 
physical  analog  of  the  limiting  behavior  of  a  thin  plate  under 
tension.  In  an  effort  to  model  the  behavior  of  the  thin  plate 
under  tension  in  two  independent  variables,  Franke  [20]  developed 
a  surface  which  is  the  analog  of  the  spline  under  tension  in  the 
same  way  that  thin  plate  splines  are  the  two  dimensional  analog 
of  cubic  splines.  The  "engineering"  approach  of  Harder  and 
Desmarais  [21]  is  followed.  The  idea  is  to  use  superposition  of 
fundamental  solutions  for  a  plate  with  tension  to  model  the 
displacement  under  point  loads.  The  fundamental  solution  (under 
appropriate  assumptions  about  the  plate  parameters;  satisfies 

42W  -  a2AW    -  6     , 
where  a    is  a  tension  parameter.   The  solution  is 


tt)    M  t  M  sKQ(as) 


r    r  t 
WQ(r)  =  (2tt)  x  /  t"1/  sKQ(as)  dsdt  +  C  , 

0   •'O 
where,   Kq  is  the  modified  Bessel  function.    The  solution  to  the 

interpolation  problem  is  obtained  by  using  this  basis  function  (C 

is   taken  to  be  zero)  in  a  global  basis  function  method  with  m=l. 

Some   experimentation  was  performed  including   linear   polynomial 

terms,  since  this  is  more  consistent  with  thin  plate  splines. 

For  tension  a   -    0 ,  the  equation  for  the  thin  plate  spline  is 

attained,   and   as  tension  gets  large,   the  equation  becomes   the 

membrane   equation,   which   has   no  finite  solution   under   point 
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loads.  Examples  are  given  which  demonstrate  that  as  tension  is 
increased  the  surface  tends  to  behave  somewhat  like  a  rubber 
sheet  under  point  loads:  too  thick  to  be  a  membrane,  but  sup- 
porting little  displacement  away  from  the  data  points. 

Construction  of  surfaces  under  tension  has  also  been  consid- 
ered by  Terzopoulous  [58].  These  surfaces  are  the  analog  of 
surfaces  constructed  by  Briggs '  method  [6],  although  the  paper 
also  addresses  the  solution  of  the  system  of  equations  by  multi- 
grid  methods.  That  appears  to  be  quite  effective,  but  will  not  be 
addressed  here.  In  addition,  as  with  Briggs'  scheme,  the  best 
situation  is  when  the  data  points  all  lie  on  a  subset  of  a 
rectangular  grid.  The  present  development  has  only  a  C  under- 
lying surface,  which  simplifies  the  calculations  for  minimization 
of  the  functional. 

The  idea  is  to  minimize  the  sum  of  a  penalty  functional 
measuring  closeness  of  fit  to  the  data  and  the  functional 


L 


p(x,y){r(x,y)(F^x+2F^y+F^y)  +  ( l-r(x, y ) ) (F^+F^ ) }dA  , 


D* 
where   D   is  the  region  of  interest  in  the  plane,   p(x,y)  is   the 

"rigidity"  of  the  plate  and  0<r(x,y)<l  is  the  "surface  tension". 
A  rectangular  grid  is  placed  over  the  region  D.  While  not  neces- 
sary in  practice,  for  simplicity  it  is  assumed  that  the  data 
points  coincide  with  grid  points.  A  nonconforming  piecewise 
quadratic  finite  element  is  assumed,  and  the  equivalent  finite 
difference  equations  for  the  above  functional  are  derived.  These 
equations  are  solved  by  multi-grid  techniques.  Local  tension  can 
be  achieved  by  allowing  the  tension  parameter  to  vary  over  the 
grid   points.    Discontinuities  in  derivative  and  value  are   also 
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considered.  As  with  Briggs'  method,  the  final  form  of  the  ap- 
proximation depends  on  the  grid  at  which  it  is  evaluated.  Unlike 
Briggs'  scheme,  which  is  based  on  finite  difference  approxima- 
tions to  the  plate  equation,  the  underlying  surface  (although  C  ) 
would  allow  evaluation  at  other  than  grid  points. 

3 . 2  Constraints 

The  construction  of  approximations  satisfying  inequality 
constraints  as  well  as  interpolation  conditions  has  been  investi- 
gated by  Villalobos  [62],  and  Dubrule  and  Kostov  [11]-  The 
former  considered  a  generalization  of  Laplacian  smoothing  splines 
(see  Section  4.1),  while  the  latter  considers  only  interpolating 
functions.  The  results  are  similar,  and  the  latter  will  be 
described  here.  The  discussion  centers  on  approximation  by  thin 
plate  splines.  The  functional  minimized  is  the  usual  thin  plate 
functional,  but  under  the  given  constraints.  As  had  been  pre- 
viously shown  by  others,  the  solution  involves  adding  basis 
functions  at  constraint  points  that  are  found  to  be  active,  and 
is  solved  as  a  quadratic  programming  problem.  This  method  ap- 
pears to  be  easy  and  effective.  The  practical  aspects  of  the 
method  are  discussed  in  the  companion  paper  [30],  where  the 
process  is  applied  using  Kriging,  which  formally  includes  thin 
plate  splines  as  a  special  case.  Journel  [28]  discusses  the 
incorporation  of  this  and  other  "soft"  information  into  the 
approximation  by  Kriging. 

A  more  general  problem  and  its  elegant  solution  is  consid- 
ered by  Utreras  [61].  Here  the  problem  is  that  of  minimizing  the 
thin  plate  functional , 
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(F2  +2F2  +F2  )  dA 
xx   xy  yy ' 

subject  to  interpolation  conditions,  and  positivity  conditions  on 
a  certain  subset  of  the  plane,  say 
F(x,y)  >  0  for  (x,y) e   K  C  R2  . 
Here  the  region  K  is  assumed  compact  and  convex.    The   existence 
and   uniqueness  of  the  solution  is  proven  using  elementary  means. 
The   characterisation  of  the  solution  of  the  problem  is  in   terms 
of   the  points  where  the  constraints   are   active.    Distribution 
theory   and   the  relationship  of  the  functional  minimized  to   the 
biharmonic   equation  are  used  to  show  that  the  solution   involves 
the    usual  terms  of  the  thin  plate  spline  approximation   plus   a 
certain   term  which  serves  to  enforce  the  positivity  of  the  func- 
tion over  the  compact  region  K. 

We  digress  to  introduce  some  notation.  Suppose  that  scat- 
tered data  is  given  with  f . > 0 .  Let  F(x,y)  be  the  solution  of  the 
constrained  problem  which  interpolates  this  data.  Define  the  set 
K*  -  {(x,y)  e  K:F(x,y)  =  0}.  K'  is  compact.  Let  B^x.y)  = 
d2log  dk,  where  d|  -  ( x-xk )2+(y-yk)2 ,  and  d2  -  x2+y2.  The 
following  theorem,  where  *  denotes  convolution,  is  then  proved. 

Theorem;    There  exist  constants  a.  and  a  positive  measure  A   with 

support  contained  in  K*  such  that 

N 
F(x,y)  =  ^ak  Bk(x,y)  +  m*BQ(x,y)  +  Px(x,y)  , 
k-1 
where  p.  is  a  polynomial  of  degree  <1 ,  and 

N 

Z-/ak  q(xk'yk}  +  m(q)  =  ° 
k=l 

for  any  polynomial  q  of  degree  <1 . 
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Convergence  is  investigated,  and  an  algorithm  is  given  for 
computing  an  approximation  to  the  positive  thin  plate  spline. 
The  algorithm  involves  finding  an  approximation  of  the  set  K*  and 
the  measure  m  such  that  the  constraint  is  satisfied  to  within 
some  tolerance.  The  approximation  is  by  points  with  atomic 
measure,  and  since  convolution  with  atomic  measures  give  point 
evaluation,  the  solution  is  approximated  by  functions  of  the  same 
type  as  in  the  Dubrule/Kostov  program.  The  crucial  difference  is 
that  here  the  constraint  set  K'  must  be  completely  discovered  as 
part  of  the  process,  rather  than  being  part  of  a  finite  subset 
specified  in  advance. 

4.0  Smoothing,  Least  Squares,  and  Subset  Selection 

In  many  applications  the  data  is  obtained  by  measurement, 
often  not  to  enough  accuracy  to  warrant  interpolation.  In  other 
applications,  such  as  oceanography  and  remote  sensing,  the  amount 
of  data  available,  even  though  subject  to  errors,  is  much  greater 
than  is  necessary  to  define  the  desired  surface  to  the  required 
accuracy.  In  these  cases  is  is  necessary  to  apply  some  smoothing 
process  or  to  perform  least  squares  approximation  with  some 
function  having  far  fewer  parameters  than  the  number  of  data 
points.  The  theory  of  smoothing  splines  in  several  variables  is 
well  developed.  For  purposes  of  contrasting  that  idea  with  those 
of  least  squares  and  subset  selection,  a  brief  discussion  of 
smoothing  splines  will  be  given. 

4.1  Laplacian  Smoothing  Splines 

Laplacian  smoothing  splines  are  a  generalization  of  univa- 
riate  smoothing   splines   to  several  variables   along   the   same 
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direction   as   thin  plate  splines  for  interpolation.    They   were 

mentioned  by  Harder  and  Desmarais  [21]  in  their  seminal  paper  for 

the   bivariate  case  under  the  physical  interpretation   of   having 

forces   applied   at  the  data  points  through  springs  with   various 

spring   constants.    This  caused  the  surface  to  tend   toward   the 

data  points,  but  the  not  necessarily  to  pass  through  them. 

The   mathematical  development  of  Laplacian  smoothing  splines 

in  the  general  case  is  given  in  Wahba  and  Wendelberger  [64] .   The 

notation  of  Section  2.1  for  points  in  s-dimensional  space  will  be 

followed   in  this  discussion.    The  functional  minimized   in   the 

case  of  smoothing  splines  is 
N 
N_1^[F(xj)-f  j]2/<7j   +  kJm^F) 
j  =  l 

where  J   is  the  pseudonorm  associated  with  interpolating  splines, 
m 

and  k  is  a  smoothing  parameter  which  governs  the  fidelity  with 
which  the  surface  fits  the  given  data.  Here  it  has  been  assumed 
that  the  errors  are  uncorrelated  and  have  standard  deviation  c. 
at  the  point  x  .  The  case  of  correlated  errors  is  addressed 
briefly  in  Wendelberger' s  thesis  [65], 

Let  a  -  (a*,  a^,  ...  ,  Of)  be  a  multi-index  for  partial 
differentiation  (denoted  by  F  ),  with  each  integer  a.  >_0  and  \d  | 
=  a^a  +    •  •  •  +as-   Then 


Jm(F)  =Z!   /  (F«)2  dx 


|at=m  •'R3 

In  this  functional  one  can  think  of  m  as  a  smoothness  parameter, 
the  approximating  functions  being  smoother  (having  higher  order 
derivatives)   the  larger  m  is.    In  order  for  the  required   func- 
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tions   to   be   in  the  appropriate  spaces,   it  is   necessary   that 

m>s/2. 

The   solution  of  the  problem  is  of  the  same  form   as   thin 

plate  splines, 
N 

F(x)  =  XX)VX)  +  qm-l(x)   ' 
j  =  l 

where   the   coefficients  a.  and  the  coefficients  ba  of  %,-a   >   a 

polynomial  of  degree  m-1,  satisfy  the  system  of  equations 
N 
Va.tB.fx1)  +  NXj^..]  +  Vbfx1)0*   =  f .  ,  i=l,2,  ...  ,N 
j=l  )a|<m 

N 
^ai(x1)0t  =  0,    |a|<m. 
1=1 

In  the  above  equations,   <5 .   is  the  Kronecker  delta, 

and  ( x)   =  X-,  X  x0   ...  x  s  . 
1     <i  s 

The  smoothing  parameter  X  must  be  specified  before  the 
smoothing  spline  can  be  computed,  and  Wahba  and  Wendelberger 
suggest  the  use  of  Generalized  Cross  Validation  (GCV)  to  select 
the  value.  Although  there  is  some  evidence  that  GCV  leads  to 
undersmoothing  (see  Seaman  and  Hutchinson  [56]),  Utreras  [60]  has 
shown  that  GCV  yields  convergence  under  reasonable  conditions. 
GCV  can  also  be  used  to  decide  on  the  smoothness  parameter,  m,  to 
be  used,  as  well.  A  limitation  of  the  scheme  is  that  it  is 
difficult  to  compute  when  there  are  more  than  200-300  data 
points,  since  the  problem  involves  the  solution  of  a  system  of 
more  than  N  linear  equations,  although  techniques  such  as  those 
of  Dyn,  described  in  the  previous  section,  could  be  useful  here. 
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4 . 2   Least  Squares  Approximation 

Least  squares  approximation  to  scattered  data  is  an  approxi- 
mation with  (presumably  many)  fewer  basis  functions  than  there 
are  data  points.  The  problem  at  hand  here,  then,  is  to  select 
the  basis  functions.  Previous  work  has  been  done  using  tensor 
product  cubic  splines,  and  a  number  of  authors  (see  ["],  [8], 
[9],  for  example)  have  considered  the  problem.  Several  computer 
programs  are  available,  and  such  methods  are  probably  desirable 
for  cases  where  data  is  somewhat  uniformly  distributed.  In  cases 
where  the  data  is  of  greatly  varying  density,  the  use  of  tensor 
products  results  in  knot  locations  on  a  grid,  and  this  may  not 
reflect  the  actual  disposition  of  data  points.  In  fact,  there 
could  be  knots  with  no  data  nearby.  While  such  problems  are  not 
insurmountable,  they  lead  to  nonuniqueness  of  the  solution,  and 
the  minimum  norm  solution  tends  to  not  be  aesthetically  pleasing. 

For  varying  density  of  data  points  it  seems  desirable  to 
have  flexibility  in  knot  placement,  and  this  leads  to  the  idea  of 
least  squares  approximation  by  thin  plate  splines.  The  MQ  ap- 
proximations also  could  be  used,  however  the  discussion  will  be 
in  terms  of  thin  plate  splines,  and  the  points  at  which  these 
basis  functions  are  centered  will  be  referred  to  as  "knot" 
points.  Treatment  of  the  knot  point  locations  as  parameters  in 
the  minimization  process  is  possible,  and  has  been  reported  upon 
by  Schmidt  [53].  The  paper  is  brief,  with  few  details  of  the 
algorithm  being  given.  The  initial  knot  configuration  was  taken 
to  be  of  tensor  product  form,  which  may  be  apparent  from  the 
final  configuration  of  knots  in  the  examples  given.  The  overall 
minimization   process   is   a  large   nonlinear   one,   and   if   the 
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problems  of  on©  dimension  carry  over  as  one  might  expect,  may  be 
complicated  in  that  knots  may  coalesce  and  the  solution  may  not 
be  unique.  In  addition  it  is  likely  true  that  the  objective 
function  may  have  many  local  minima,  so  an  algorithm  to  search 
for  a  better  local  minimum,  or  to  avoid  a  poor  local  minimum 
would  be  necessary. 

The  problem  with  treating  the  knot  locations  as  parameters 
presently  seems  to  be  intractable  mathematically,  and  for  many 
knots  is  computationally  expensive,  with  results  obtained  being 
of  questionable  quality.  A  different  point  of  view  on  the  prob- 
lem is  considered  in  McMahon  [36],  and  a  summary  of  his  approach 
and  results  will  be  given. 

The  main  problem  in  the  process  is  the  selection  of  knots, 
for  once  these  have  been  decided  upon,  what  remains  is  to  solve 
the  equations  in  the  least  squares  sense,  in  principle  an  unin- 
teresting task.  If  the  selection  of  knot  locations  is  to  be 
decoupled  from  the  least  squares  process,  some  assumption  must  be 
made  order  to  have  an  algorithm  for  selection  of  the  knots.  The 
assumption  in  this  case  is  that  the  independent  variable  data 
reflects  something  about  the  behavior  of  the  dependent  variable. 
For  example,  perhaps  the  density  of  data  point  locations  is 
dependent  on  the  curvature  of  the  surface,  or  more  broadly,  if 
the  function  is  changing  behavior  rapidly  the  density  of  data 
points  is  great,  whereas  low  density  indicates  slowly  changing 
behavior.  This  assumption  is  not  universally  satisfied  in  prac- 
tice, for  some  data  is  taken  based  on  accessibility  (along  roads 
in  rugged  areas,  for  example)  or  other  nonbehavioral  criteria. 

The  assumption  of  data  density  indicating  local  behavior   of 
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the  surface  has  a  companion  assumption  that  each  data  point  is  in 
some  sense  equally  important  in  defining  the  function,  regardless 
of  whether  it  has  a  nearby  neighbor  or  not.  This  leads  to  the 
idea  of  "equal  representation"  for  each  data  point  by  a  knot 
point.  This  means  that  each  data  point  should  be  "close"  to  a 
knot  point,  and  that  each  knot  point  should  "represent"  about  the 
same  number  of  data  points.  These  two  ideas  are  crucial  to  the 
knot  selection  algorithm.  First,  it  is  desirable  to  minimize  the 
sum  of  the  distances  squared  from  each  data  point  to  the  nearest 
knot  point.  Let  the  knot  point  locations  be  given  by  (St.,f.), 
j-1 ,     2,     . . . ,    K.   Then  the  function  to  be  minimized  is 


N 


GN2  =  ^minnx^-a^+Cy^-y:,)2; 
k=l  j 


This  process  has  a  default  Dirichlet  tesselation  with  respect  to 
the  knot  points,  with  each  data  point  "belonging"  to  some  knot 
point  by  virtue  of  the  Dirichlet  tile  to  which  it  belongs.  When 
a  data  point  lies  on  a  tile  boundary,  some  determination  of  which 
it  belongs  to,  or  whether  it  is  shared  must  be  assumed.  Local 
minima  for  the  above  function  are  easily  characterized:  At  a 
local  minimum,  each  knot  point  (x.,y.)  is  the  centroid  of  the 
data  points  in  its  tile.  This  leads  to  a  nice  algorithm  for 
iteration  to  a  local  minimum,  by  repeated  computation  of  the 
centroid  of  the  data  points  in  the  Dirichlet  tiles. 

This  algorithm  works  very  well  for  obtaining  a  local  minimum 

9 

of   the  function  GN  .    Unfortunately  the  function  is   rife   with 

local  minima,  so  that  finding  a  desirable  one  depends  on  the 
proper  initial  guess.  Here  the  second  idea  regarding  "equal 
representation"   comes   into   play.    Since  each   data   point   is 
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assumed  equally  important,  it  is  reasonable  to  expect  that  the 
Dirichlet  tile  for  each  knot  should  contain  about  the  same  number 
of  data  points.  This  leads  to  a  heuristic  for  a  new  initial 
guess   at   knot  locations  once  a  local  minimum  for  GN    has   been 

found.   It  is  still  considered  desirable  to  be  at  a  local  minimum 

2  2 

of   GN  ,   but   the   idea  of  finding  a  global  minimum   of   GN    is 

abandoned   for   the  "equal  representation"  idea.    Once   a   local 

minimum  has  been  found,   a  new  measure  of  "goodness"  is  computed. 

Let   N.  be  the  number  of  points  in  the  Dirichlet  tile  for  the  j 

knot  point  (£.,y.).   Let 


K 


-E 


(N^N/K)2  , 


j  =  l 

a   measure   of  the  "equal  representation"  for  a   particular   knot 

configuration.   It  was  then  attempted  to  determine  knot  locations 

2 
which  achieve  a  local  minimum  of  GN   and  which  also  yield  a  small 

value,  hopefully  a  minimum,  of  D. 

The  general  idea  of  the  algorithm  is  to  "nudge"  the  knot 
locations  from  a  local  minimum  of  GN  toward  a  configuration  with 
smaller  D  value.  The  rationale  for  this  is  to  attempt  to  move 
the  tile  boundaries  across  data  points,  yielding  a  more  equitable 
distribution.  The  search  for  the  minimum  of  D  could  turn  out  to 
be  rather  extensive,  and  for  a  large  number  of  points  with  a 
moderately  large  number  of  knot  points,  the  computational  effort 
can  be  excessive.  End  results  are  still  somewhat  dependent  on 
the  initial  guess,  although  the  results  generally  tend  to  look 
quite  reasonable.  The  program  incorporates  the  option  of  inter- 
nally generated  ( quasi-gridded)  or  user  input  initial  guesses. 
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Examples  are  given  which  illustrate  rather  well  the  ability 
of  the  scheme  to  select  knot  locations  which  reflect  the  under- 
lying density  of  the  data.  Actual  surface  fitting  and  comparison 
with  two  other  methods,  the  Laplacian  smoothing  splines  of  Wahba 
and  Wendelberger ,  and  the  tensor  product  bicubic  Hermite  method 
due  to  Foley  [18],  are  reported  upon.  One  example  illustrates 
the  failure  of  the  fitting  scheme  to  accurately  model  the  surface 
when  the  data  density  does  not  reflect  the  underlying  behavior  of 
the  surface,  as  was  assumed. 

9 

4 . 3   Subset  Selection 

To  my  knowledge,  the  first  investigation  into  the  use  of  a 
subset  of  scattered  data  points  upon  which  to  base  an  interpolant 
to  be  used  as  an  approximation  for  the  entire  data  set  was  by 
Pickrell  [47] .  The  application  guided  many  of  the  ideas  involved 
in  the  thesis  and  they  are  not  universally  applicable. 
Pickrell 's  goal  was  to  model  underwater  terrain  from  a  large 
number  of  measured  depths,  which  were  quite  accurate,  the  most 
important  being  those  which  were  on  ridges  or  other  shallow 
areas,  since  the  resulting  approximation  could  be  used  to  gener- 
ate charts  for  navigation. 

The  idea  was  to  select  as  small  a  subset  of  points  as  pos- 
sible such  that  the  interpolant  for  these  points  yielded  a  sur- 
face that  satisfied  a  certain  error  tolerance  at  all  the  other 
data  points.  An  iterative  scheme  was  used.  Beginning  with  an 
initial  guess,  either  taken  "uniformly"  spaced  over  the  region  of 
interest,  or  by  perusal  of  the  data  for  "critical"  points  such  as 
high  and  low  values  or  areas  of  sharp  gradients.   Call  the  selec- 
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ted  points  the  "model"  points.  The  interpolant  (the  MQ  method 
was  used;  parameter  r=0  generally  performed  best)  was  constructed 
and  the  deviation  at  other  data  points  computed.  Based  on  this 
information,  other  points  at  which  the  deviation  was  large,  or 
one  from  a  group  of  nearby  points  where  deviations  of  the  same 
sign  occurred,  would  be  included  in  the  set  of  model  points. 
Points  with  small  coefficients  might  be  eliminated  from  the  set 
of  model  points,  as  well.  The  investigation  did  not  yield  an 
algorithm  suitable  for  approximation  with  no  intervention  by  the 
user,  since  this  was  performed  interactively  by  Pickrell.  In  the 
examples  reported,  the  required  number  of  model  points  was  gen- 
erally on  the  order  of  10-20%  of  the  total  number  of  points. 
However,  only  limited  consideration  was  given  to  the  handling  of 
very  large  data  sets  which  would  require  local  application  of  the 
ideas  with  a  scheme  for  joining  the  pieces  together. 

These  deficiencies  and  other  matters  were  addressed  by 
Schiro  and  Williams  [52].  In  addition  to  automating  the  model 
point  selection  process,  he  subdivided  the  region  into  kernel 
groups  based  on  a  measure  of  homogeneity  of  the  data.  As  a 
starting  point  the  region  of  interest  was  divided  into  a  rectang- 
ular grid  of  cells  of  size  equal  to  the  minimum  to  be  considered. 
The  mean  and  standard  deviation  of  the  function  values  (again, 
ocean  depth  data)  were  computed  for  each  cell.  Then,  the  cells 
were  considered  in  turn  as  a  base  cell,  and  contiguous  cells 
(with  larger  coordinates)  having  mean  values  within  one  standard 
deviation  of  the  mean  for  the  base  cell  were  combined  to  form  a 
larger  rectangular  cell,  when  possible.  This  larger  group  of 
cells  is  called  a  kernel  group.    Thus,   each  kernel  group  in  the 
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final  subdivision  was  made  up  of  one  or  more  contiguous  cells  in 
the  original  subdivision,  all  of  which  have  similar  mean  behavior 

The  selection  of  model  points  for  each  kernel  group  is  done 
without  user  intervention.  The  initial  set  is  chosen  based  on 
deviation  from  the  mean  value.  (Of  course,  if  all  points  are 
within  the  tolerance  of  the  mean  value,  the  mean  value  is  used  as 
the  approximation  for  the  kernel  group. )  The  idea  is  similar  to 
Pickrell's:  points  with  large  deviation  are  selected  to  be  added 
to  the  set  of  model  points,  with  the  proviso  that  points  closer 
than  a  certain  distance  cannot  be  selected  on  the  same  iteration. 
This  avoids  the  problem  of  adding  several  points  very  close 
together  on  the  same  iteration.  The  process  is  continued  until 
the  deviations  are  below  a  specified  tolerance. 

The  overall  surface  is  made  continuous  by  blending  adjacent 
surfaces  when  within  a  certain  distance  from  boundaries  of  the 
kernel  groups.  This  is  done  via  Hermite  cubics  to  obtain  a 
smooth  transition.  The  MQ  method  with  r=0  was  used  to  fit  the 
differences,  data  minus  mean  value,  and  this  was  observed  to  have 
a  beneficial  effect  in  terms  of  the  magnitude  of  the  coefficients 
in  the  representation. 

Another  approach  to  subset  selection  was  taken  by  Bozzini , 
deTisi,  and  Lenarduzzi  [5].  Here  an  attempt  was  made  to  deter- 
mine a  subset  of  points  to  be  used  to  compute  the  approximation 
by  doing  local  computations  to  determine  whether  a  particular 
point  has  a  significant  influence  on  defining  the  surface.  A 
description  of  the  ideas  will  be  given.  The  first  step  is  to 
partition  the  region  of  interest  into  (probably  overlapping) 
regions   of   given  shape  (e.g.,   circular  disks)   by   taking   the 
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region  R.  as  large  as  possible  so  that 


77^-t/  lf(x)-f. 
a(RiVR.       i 


,|  dx  <  K. ,  for  some  constant  K. 
area  Ki  )  / •»       j-        1  1 


The  value  of  K.  specifies  a  kind  of  homogeneity  of  f(x)  over  R. . 
Then,  for  a  weighting  function  c$ .  (x),  let 

d.    =  /  jf(x)-f.|  <*(x)  dx  //  4>(x)  dx 

This  value  serves  as  measure  of  the  level  of  homogeneity  in  the 
region  R .  ,  weighted  by  4> .  Then  a  measure  of  the  behavior  of  f(x) 
at  point  i  relative  to  point  j  (point  j  also  in  R. )  is  given  by 

These  values  are  averaged  over  the  points  in  R.  to  obtain  s. . 
The  value  of  3.  is  a  measure  of  the  what  the  authors  call  the 
"importance"  of  the  point  i  in  the  data  set,  and  hence  to  the 
definition  of  the  approximating  surface.  The  integrals  in  the 
definition  are  approximated  with  the  obvious  quadratures  in  the 
application.  The  point  with  the  largest  corresponding  value  of 
s.  is  chosen.  The  process  can  then  be  repeated  from  the  computa- 
tion of  the  s.  to  obtain  more  points,  or  one  can  choose  the 
subset  with  the  largest  s.  values  from  the  initial  calculation. 
According  to  the  authors,  the  method  is  not  to  sensitive  to 
deletion  of  a  point  from  the  set.  Some  examples  are  given  to 
illustrate  approximation  of  surfaces  from  subsets  of  a  given  set 
of  points. 
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